(*
** some code for testing various functions in gsl_statistics
*)

staload UN = "prelude/SATS/unsafe.sats"


staload "libc/SATS/time.sats"

staload "contrib/gsl/SATS/gsl_rng.sats"
staload "contrib/gsl/SATS/gsl_randist.sats"
staload "contrib/gsl/SATS/gsl_statistics.sats"


staload "prelude/SATS/array.sats"
staload "prelude/DATS/array.dats"


extern
fun get_numbers{n:nat}{l:agz}
  (r: !gsl_prng l,cnt: int n) : list(double,n)

implement get_numbers{n} (r,cnt) = 
  if cnt = 0 then
    list_nil
  else 
    let
      (* Get a zero mean guassian with standard dev. 5*)
      val standard = gsl_ran_gaussian(r,5.0)
      (* Give the guassian a mean of 10.*)
      val result = standard+10.0
    in
      list_cons(result,get_numbers(r,cnt-1))
    end

#define N 100000

implement main () = let
  (* Setup RNG *)
  val _ = gsl_rng_env_setup ()
  val my_rng = gsl_rng_alloc_exn (gsl_rng_default)
  val t1 = time_get ()
  val () = gsl_rng_set (my_rng, $UN.cast2ulint (t1))
  
  (* Setup an array of values. *)
  val (pfr,pf | results) = array_ptr_alloc_tsz{double}(N,sizeof<double>)
  val () = array_ptr_initialize_lst<double>(!results,get_numbers(my_rng,N))
  
  (* Collect Some Information *)
  val mean = gsl_stats_mean(!results, 1, N)
  val variance = gsl_stats_variance(!results,1, N)
  val sd = gsl_stats_sd(!results,1,N)
  val skew = gsl_stats_skew(!results,1,N)
  val kurtosis = gsl_stats_kurtosis(!results, 1,N)
  
  (* Display Values *)
  val () = printf("Mean is: %lf\n",@(mean))
  val () = printf("Variance is %lf\n",@(variance))
  val () = printf("Standard Deviation is: %lf\n",@(sd))
  val () = printf("Skew is: %lf\n",@(skew))
  val () = printf("Kurtosis is: %lf\n",@(kurtosis))
  val () = array_ptr_free(pfr,pf | results)
  val () = gsl_rng_free(my_rng)
  in end